optedr: un paquete de diseño óptimo de experimentos


Carlos de la Calle Arroyo

https://github.com/Kezrael/R-Hispano-2022-optedr

cdelacallea@unav.es

Institute of Data Science and Artificial Intelligence

24/11/2022

Diseño Óptimo de Experimentos


Modelo de regresión

\[y(x)=\eta(x;\theta)+\varepsilon\quad, \]

con \(\theta\) los parámetros desconocidos del modelo, de interés para el experimentador.

Diseño Experimental

  • Diseño exacto: una colección de puntos, con sus réplicas, en los que tomar medidas.
  • Diseño aproximado: un diseño, \(\xi\), visto como una medida de probabilidad.

FIM

La Matriz de Información de Fisher se define como

\[ M(\xi,\theta)=\sum_{x\in\mathcal{X}}f(x)f^t(x)\xi(x), \]

con \(f(x)=\partial\eta(x, \theta)/\partial\theta\)

Criterios de optimalidad

  • D-optimalidad: \(\Phi_D=|M^{-1/k}|\).
  • \(D_s\)-optimalidad: \(\Phi_{D_S}=\left \{\frac{|M|}{|M_{22}|}\right\}^{-1/s}\).
  • A-optimalidad: \(\Phi_A=\text{Tr} (M^{-1})\).
  • I-optimalidad: \(\Phi_I=\text{Tr} (AM^{-1})\), con \(A=\int_{\mathcal{R}}\mu(x)f(x)f(x)^Tdx\).

Teorema General de Equivalencia

Un diseño \(\xi^\star\) es \(\phi\)-óptimo si y sólo si,

\[Tr \left\{\nabla \phi(\xi^\star)[M(\xi^\star, \theta)-f(x)f^t(x)]\right\}\leq 0, x\in\mathcal{X}.\]

La igualdad se cumple en los puntos de soporte del diseño.

Eficiencia

La eficiencia de un diseño \(\xi\), para un cierto cirterio \(\Phi\), se define como

\[ eff_\Phi(\xi) = \frac{\Phi(\xi)}{\Phi(\xi^\star)} \]

Cálculo de diseños óptimos


Algoritmo Cocktail

  1. Empezar con un diseño no singular.
  2. Añadir \(x\) tal que \(\max_{x\in\chi} f^t(x)\nabla \phi(\xi)f(x)\).
  3. Optimizar los pesos a través del algoritmo multiplicativo.
  4. Aplicar heurísticos.
  5. Comprobar optimalidad, si no, repetir 2-5.

Cálculo de diseños óptimos: función opt_des

resArr.D <- opt_des("D-Optimality",
  model = y ~ a*exp(-b/x),
  parameters = c("a", "b"),
  par_values = c(1, 1500),
  design_space = c(212, 422))

Output

kable(resArr.D$optdes)
Point Weight
329.2966 0.5000068
422.0000 0.4999932
resArr.D$sens

Función de sensibilidad del diseño

resArr.D$convergence

Convergencia del algoritmo

Distribuciones de probabilidad

Se puede seleccionar con el parámetro distribution

  • ‘Homoscedasticity’.
  • ‘Gamma’, que puede ser usada también para una exponencial o una normal heterocedástica con error relativo constante.
  • ‘Poisson’.
  • ‘Logistic’.

también se puede dar la función directamente con el parámetro weight_fun.

Diseño para la distribución de Poisson.

poissionDopt <- opt_des("D-Optimality",
   y ~ t0 + t1 * x + t2 * x^2, 
   c("t0","t1","t2"),
   c(0.95, -1, -1), 
   c(0, 0.5799), 
   distribution = "Poisson")
Point Weight
0.0000000 0.3333275
0.2900226 0.3333274
0.5799000 0.3333451

Heurísticos del algoritmo

FP2_01 <- opt_des("I-Optimality",
  y ~ a0 + a1*log(x) + a2*log(x)^2,
  c("a0", "a1", "a2"),
  c(1, 1, 1),
  c(0.0001, 1),
  reg_int = c(0.0001, 1))
Point Weight
0.000387 0.0441516
0.166750 0.4951323
1.000000 0.4607161

que tiene una eficiencia de \(81.83\%\).

Heurísticos del algoritmo

FP2_02 <- opt_des("I-Optimality",
  y ~ a0 + a1*log(x) + a2*log(x)^2,
  c("a0", "a1", "a2"),
  c(1, 1, 1),
  c(0.0001, 1),
  reg_int = c(0.0001, 1),
  join_thresh = 0.005)
Point Weight
0.0001000 0.0366361
0.0335051 0.3767898
1.0000000 0.5865741

con una eficiencia mayor a \(99.99\%\).

Evaluando diseños: la clase optdes


Eficiencia de un diseño

Un diseño cualquiera

design <- data.frame("Point" = c(220, 240, 400),
  "Weight" = c(1/3, 1/3, 1/3))

design_efficiency(resArr.D, design)
[1] 0.3063763

O los dos diseños calculados antes

design_efficiency(FP2_02, FP2_01$optdes)
[1] 0.8183762

optdes

print(resArr.D)
     Point    Weight
1 329.2966 0.5000068
2 422.0000 0.4999932
summary(resArr.D)
Model: 
y ~ a * exp(-b/x)
and weight function: 
NULL
Optimal design for  D-Optimality :
     Point    Weight
1 329.2966 0.5000068
2 422.0000 0.4999932

 Minimum efficiency (Atwood):  99.9986396401789%
 Criterion value:  9972806

Clase optdes

plot(resArr.D)

Plot del diseño óptimo

Diseños aumentados


Región de puntos candidatos

get_augment_region("D-Optimality",
               init_design = resArr.D$optdes, 
               alpha = 0.3,
               model = y ~ a * exp(-b/x),
               parameters = c("a", "b"),
               par_values = c(1, 1500),
               design_space = c(212, 422),
               calc_optimal_design = F)

Región de puntos candidatos

Diseños aumentados

aug_arr <- augment_design(
               "D-Optimality",
               init_design = resArr.D$optdes, 
               alpha = 0.3,
               model = y ~ a * exp(-b/x),
               parameters = c("a", "b"),
               par_values = c(1, 1500),
               design_space = c(212, 422),
               calc_optimal_design = F)

Eficiencia de \(80\%\)

Añadiendo los puntos 260 y 380 con peso 0.15

aug_arr
     Point    Weight
1 329.2966 0.3500048
2 422.0000 0.3499952
3 260.0000 0.1500000
4 380.0000 0.1500000

Redondeo eficiente


Algoritmo de redondeo

  1. Empezar con un diseños aproximado con \(k\) puntos y pesos \(\xi(x_i)\)
  2. Calcular \(\nu = n - k/2\)
  3. Asignar \(z_i \gets \nu\xi(x_i)\)
  4. Redondear \(z_i\) al entero más cercano \(r_i = \lceil {z} \rceil\)
  5. Si \(\sum_i r_i = n\) entonces devolver \(r = \{r_1, \dots, r_n\}\) si no incrementar \(r_i\) con \(i = argmax_i (n \xi(x_i) - r_i)\) o disminuir \(r_i\) con \(i = argmin_i (n \xi(x_i) - r_i)\) hasta \(\sum_i r_i = n\)

Redondeo

design_test <- data.frame("Point" = seq(1, 5, length.out = 7),
  "Weight" = c(0.1, 0.0001, 0.2, 0.134, 0.073, 0.2111, 0.2818))

efficient_round(design_test, 21)
Point Weight
1.000000 2
1.666667 1
2.333333 4
3.000000 3
3.666667 2
4.333333 4
5.000000 5

Comprobar eficiencia

exact_design <- efficient_round(resArr.D$optdes, 7)
aprox_design <- exact_design
aprox_design$Weight <-
  aprox_design$Weight / sum(aprox_design$Weight)

design_efficiency(resArr.D, aprox_design)
[1] 0.9897433

Shiny demos


Shiny Antoine

shiny_optimal()

Shiny de diseños óptimos para la Ecuación de Antoine

Shiny Addpoints

shiny_augment()

Shiny de aumentar diseños

Resumen y desarrollo futuro


Resumen

  • Facilidad de uso
  • Aumentar diseños para diferentes modelos
  • Redondeo eficiente para poder implementar los diseños

Futuras mejoras

  • Modelos con varias variables explicativas
  • Interfaz gráfica con funcionalidad completa

El paquete en github y CRAN

  • https://github.com/Kezrael/optedr
  • https://cran.r-project.org/web/packages/optedr/index.html

¡Gracias por vuestra atención!